/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "uhm.hxx"

#define UHM_ERROR_TOL 1.0e-5

int main (int argc, char **argv)
{
  // input check
  if (argc != 4) {
    printf("Try : %s [decomposition][delta_dim][input_file]\n\n", argv[0]);
    printf(" ==== decomposition method ====\n\n");
    printf(" ** Symmetric Factorization **\n");
    printf("    UHM_CHOL      ( 1),   UHM_CHOL_PIV ( 2),   UHM_CHOL_INCPIV ( 3) \n");
    printf("    UHM_LDL_NOPIV ( 4),   UHM_LDL_PIV  ( 5),   UHM_LDL_INCPIV  ( 6) \n");
    printf("\n");
    printf(" ** General Structure Symmetric Factorization **\n");
    printf("    UHM_LU_NOPIV  ( 7),   UHM_LU_PIV   ( 8),   UHM_LU_INCPIV   ( 9) \n");
    printf("    UHM_QR        (10) \n");
    return 0;
  }

  FLA_Init();

  uhm::Mesh m;

  int decomposition, delta_dim;
  char *filename;

  decomposition = atoi( (argv[1]) );
  delta_dim     = atoi( (argv[2]) );
  filename      = argv[3];

  printf( "BEGIN : Import mesh from file : %s\n", filename );
  m = new uhm::Mesh_;
  m->import_file( filename );
  printf( "END   : Import \n" );

  uhm::set_num_threads(8);

  printf( "BEGIN : Build tree \n" );
  uhm::build_tree(m);
  printf( "END   : Build tree\n" );

  printf( "New mesh info < n_nodes %d, n_elements %d >\n", 
	  m->get_n_nodes(), m->get_n_elements() );

  m->lock();

  int datatype = UHM_REAL, n_rhs = 1, is_schur = false;

  uhm::set_hier_blocksize(256);

  if (decomposition < UHM_LU_NOPIV)
    m->set_symmetry(UHM_SYMMETRY);
  else
    m->set_symmetry(UHM_UNSYMMETRY);
  printf( "BEGIN : Create matrix without buffer \n" );
  m->create_matrix_without_buffer( datatype, n_rhs );
  printf( "END   : Create matrix without buffer \n" );

  // ----------------------------------------------------------------
  unsigned int n_dof     = m->get_n_dof();
  double f_decompose, f_solve, m_estimate;
  m->estimate_cost( decomposition, datatype, n_rhs,
                    f_decompose, f_solve, m_estimate );

  std::vector<int> x, dist_uhm;
  std::vector<double> dist_flop, dist_buffer;
  m->get_matrix_stat( decomposition, datatype,
                      delta_dim, x, dist_uhm, dist_flop, dist_buffer );

  printf("==== Histogram of UHM =====\n");
  printf("dim   freq   flop   buffer \n");
  for (int i=0;i<x.size();++i) 
    printf("%d  %d  %E %E \n", x[i], dist_uhm[i], dist_flop[i], dist_buffer[i]);


  printf("==== Report =====\n");
  printf("Decomposition          = %d\n", decomposition);
  printf("NDOF                   = %d\n", n_dof);
  printf("Buffer estimate (MB)   = %6.0lf\n", m_estimate/1.0e6);
  printf("FLOP decom (GFLOP)     = %6.3lf\n", f_decompose/1.0e9);
  printf("FLOP solve (GFLOP)     = %6.3lf\n", f_solve/1.0e9);
  printf("--------------------------\n");

  m->unlock();

  delete m;

  FLA_Finalize();
  return 0;
}


